library(foreign)

PRD <- function(x,P){mean(x)*mean(P)/mean(P*x)}

astats = function(av,price){
  ratio = av/price
  rmean = mean(ratio)
  rmedian = median(ratio)
  wgtmean = mean(price*ratio)/mean(price)
  cod = 100*mean(abs(ratio - rmedian)/rmedian)
  prd = rmean/wgtmean
  out = c(rmean, rmedian, wgtmean, cod, prd)
  return(out)
}

Suits <- function(av,price,nboot=0){
  data <- data.frame(av,price)
  o <- order(data$price)
  data <- data[o,]
  cumav <- 100*cumsum(data$av/1000)/sum(data$av/1000)
  cumprice <- 100*cumsum(data$price/1000)/sum(data$price/1000)
  n = length(cumav)
  dav <- cumav[2:n] + cumav[seq(1,n-1)]
  dav=c(cumav[1],dav)
  dprice <- cumprice[2:n] - cumprice[seq(1,n-1)]
  dprice=c(cumprice[1],dprice)
  L = sum(.5*dav*dprice)
  S =  1 - L/5000
  out = S

  if (nboot>0){
    outse <- array(0,dim=nboot)
    for (i in 1:nboot){
      obs <- sample(1:n,n,replace=T)
      newdata <- data[obs,]
      o <- order(newdata$price)
      newdata <- newdata[o,]
      cumav <- 100*cumsum(newdata$av)/sum(newdata$av)
      cumprice <- 100*cumsum(newdata$price)/sum(newdata$price)
      dav <- cumav[2:n] + cumav[seq(1,n-1)]
      dav=c(cumav[1],dav)
      dprice <- cumprice[2:n] - cumprice[seq(1,n-1)]
      dprice=c(cumprice[1],dprice)
      L = sum(.5*dav*dprice)
      S =  1 - L/5000
      outse[i] = S
    }
    se = sd(outse)
    out <- c(out,se)
  }
  return(out)
}

ExemptSuits <- function(E,m,av,price,g=1){
  taxbase <- pmax(0,m*av - E)
  out <- abs(Suits(taxbase,price))
  return(out)
}

ExemptPRD <- function(E,m,av,price,g=1){
  taxbase <- pmax(0,m*av - E)
  out <- abs(PRD(taxbase/price,price)-g)
  return(out)
}


rdata <- read.dta("./CoreLogic/ratio/ratio.dta")
table(rdata$year,rdata$avyear)
sampvar <- abs(rdata$avyear-rdata$year)<=2
rdata <- rdata[sampvar,]
rdata$fipscity <- paste(rdata$fips,rdata$city)
fipscity <- unique(rdata$fipscity)
n = length(fipscity)
n
nobs <- array(1,dim=nrow(rdata))
fit <- aggregate(nobs,by=list(rdata$fipscity),sum)
names(fit) <- c("fipscity","nobs")
fit100 <- fit[fit$nobs<100,]
nrow(fit100)
rdata <- rdata[!rdata$fipscity%in%fit100$fipscity,]
dim(rdata)
fipscity <- unique(rdata$fipscity)
n = length(fipscity)
n
table(rdata$year,rdata$avyear)



outdata <- NULL
for (i in 1:n){
  if (i%in%seq(0,n,100)){cat(i,"\n")}
  sampvar <- rdata$fipscity==fipscity[i]
  av <- rdata$av[sampvar]
  price <- rdata$price[sampvar]
  ratio <- av/price
  suits = round(Suits(av,price),3)
  fips = rdata$fips[sampvar][1]
  city = rdata$city[sampvar][1]
  avstats <- round(astats(av,price),3)
  mprice = round(median(price))
  
  out <- c(fips,city,length(av), avstats,mprice,suits)
  prd = avstats[5]

  m = 1/median(ratio)
  
  if (suits<0){
    for (i in seq(10000,5000000,10000)){
      E <- round(optimize(ExemptSuits,c(i-10000,i),m,av,price)$minimum)
      if (E!=i){break}
    }
  }  else {E = 0}
  out <- c(out,E)
  if (prd>1.0){
    for (i in seq(10000,5000000,10000)){
      E <- round(optimize(ExemptPRD,c(i-10000,i),m,av,price)$minimum)
      if (E!=i){break}
    }
  }  else {E = 0}
  out <- c(out,E)

  outdata <- rbind(outdata,out)
}
colnames(outdata) <- c("fips","city","nobs","Mean","Median","wmean","COD","PRD","mprice","Suits","Exempt1","Exempt2")
rownames(outdata) <- 1:nrow(outdata)
outdata <- data.frame(outdata)
for (j in 4:12){
  outdata[,j] <- as.numeric(outdata[,j])
}

write.dta(outdata,file="./research/assess/Exempt/city.dta",version=12)




